The new evolutionary code Pegase.3 is used for the spectral synthesis of continuous optical-Spitzer-Herschel-submm spectral energy distributions (SEDs) of two z=3.8 radio galaxies 4C 41.17 and TN J2007-1316. Both galaxies are selected from the Herschel Radio Galaxy Evolution Project for their faint AGN contribution and because they show evidence of a large stellar contribution to their bolometric luminosity. Pegase.3 coherently models the reprocessing of the stellar luminosity to dust emission, allowing to build UV to IR-submm libraries by types. These templates allow to predict SEDs at redshifts z in the observer's frame by using z=0 templates corrected for cosmology and evolution. Our best fits are a sum of two distinct evolving stellar populations: i) a very massive starburst observed ~30 Myrs after formation from optical and far-IR peaks ii) a significantly older (~1Gyr old) massive early-type population revealed from mid-IR Spitzer data. The AGN emission fits a roughly thermal model. These results confirm that many of the stellar populations in high-z radio galaxies were already formed by massive (a few thousands of billion solar masses) starbursts at z=4, confirming the previous K-z diagram interpretation. Gas-rich mergers and/or jet-cloud interactions are favored for triggering the intense star formation, possibly in relation with massive super black holes of these radio galaxies. These multiple stellar populations may be a generic feature of the luminous infrared radio galaxy population as a whole and strongly constrain mass accumulation and galaxy formation models.